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Abstract 

The interactions between charge and orbitaUy ordered d-electrons are important in many tran- 
sition metal oxides. We propose an effective energy model for such interactions, parameterized 
with DFT+U calculations, so that energy contributions of both electronic and lattice origin can 
Q be simultaneously accounted for. The model is applied to the low-temperature phase of magnetite, 

for which we propose a new ground state structure. The effective interactions on the B-lattice 
of Fe304 can be interpreted in terms of electrostatics and short-range Kugel-Khomskii exchange 
coupling. The frustration between optimal charge and orbital orderings leads to a complex energy 
landscape whereby the supercell for the charge ordering, orbital ordering and ionic displacements 
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PACS numbers: 71.10.Li, 71.28.+d, 75.30.Et 



I. INTRODUCTION 

The physical properties of transition-metal oxides (TMO) are determined by two degrees 
of freedom that are often closely related: the quantum state of interacting electrons, and the 
ionic position and/or motion. By itself the description of correlated (i-electrons, if bound 
to or nearly localized on specific transition metal sites, is already complex. The magnetic, 
charge (electron or hole) and orbital (when the d shell of the transition metal ion is partially 
filled) degrees of freedom of localized d-electrons are coupled via electrostatic, direct and/or 
super-exchange interactions (for a review, see e.g. refs. HHSj). The interplay between these 
degrees of freedom and their possible ordering play an important role in understanding such 
phenomena as metal-insulator transitions, high-temperature superconductivity and colossal 
magnetoresistance. On the other hand, the ionic displacements may mediate Jahn- Teller 
interactions of degenerate orbitals and induce orbital ordering (00) (see e.g. ref. H]). The 
energetic effects of both mechanisms are often of the same order of magnitude (~ 10 — 
10^ meV). It is therefore of theoretical and practical interest to investigate the combined 
electronic and lattice effects on the electronic ground state. 

First-principles calculations based on the density functional theory (DFT) provide a nat- 
ural way to incorporate both electronic and ionic degrees of freedom in "real-world" mate- 
rials. Since such methods provide direct information only about the energy of the system 
as a whole, the microscopic relationship between the involved degrees of freedom has to be 
extracted indirectly. 

Recently, we showed how the effective interactions of localized d-electrons (minority 
spin d-state tLe^ of high-spin Fe^"*") and holes (high-spin Fe^"*") in the mixed-valence oxide 
Li-EFeP04 could be extracted from first-principles calculations^. Since two types of charge 
carriers, Li-ions and d-electrons, coexist in this material, our energy model includes the ionic, 
ion-electron interactions as well. The inter-site coupling parameters could be obtained with 
the cluster expansion approach^ whereby an Ising-like Hamiltonian in electron occupation 
variables is fitted to DFT-I-U— total energy calculations of different charge ordering (CO) 
and ionic ordering patterns. It was found that the effective electron interactions exhib- 
ited strong electrostatic character at short range and lattice influence at long range. The 
accuracy of this approach was supported by the good agreement between the computecP 
and experimentaP^ temperature-composition phase diagram. Since the t2g degeneracy in 



Li-EFeP04 is lifted by the irregular FeOe octahedra, we did not explicitly consider the orbital 
degree of freedom for Fe^"*" in that work. 

In a system with degenerate localized d or f states, multiple self-consistent Kohn-Sham 
solutions may appear, a reflection of the existence of orbital ordering. Since / electrons are 
well localized, calculations are often trapped in local minima, making it difficult to finding 
the ground state (for very recent discussions, see vei!^. In comparison, the strong crystal 
field effects in d-systems make the calculations relatively straightforward. The computational 
methods and settings are presented in section |Tl| the obtained results and relevant discussions 



in section III, and finally the conclusions and outlook for future work in section IV 



II. COMPUTATIONAL METHOD 



A. Energy model 



In this paper we explore an approach to derive orbital physics from DFT total energies. 
A model is introduced for first-principles determination of the effective interactions of charge 
and orbital oideied (COO) electrons, and applied to Fe304. The model includes electrostat- 
ics, lattice distortion etc in a consistent manner. Consider a general energy expression: 

E[e\ = Eo + Ei{ei) + Eij{ti, e^) + Eijk{ei, e^, e^) H (1) 

where e^ represents the electronic state (hole and/or orbital) on site i and e'is the system's 
configuration of states. Summation over repeating indices are implied. The point term 
Ei describes the electron chemical potential and splitting of the orbitals. Besides the pair 
interaction matrices Eij, one includes in general higher order contributions, e.g. Eij^. In 
practice this model may be too general to use. Preempting the finding that quantum ef- 
fects that distinguish the different orbitals vanish at long distance, it is more convenient to 
separate orbital-independent interactions Ec from orbital-dependent Eq, attributed to the 
charge and the orbital degrees of freedom, respectively. The former can be described by a 
binary (electron/hole) cluster expansion modeP'^, which captures both short and long-range 
effects^. Therefore we rewrite eq. [l] as 

E[e] = E,[e] + E,[e] (2) 

Ec[e\ = Ec[e] = J© + Mi + Jijiiij H 



where J's are effective cluster interactions (ECI) of a cluster expansion^, e = — 1 when e is the 
hole state and +1 otherwise, and V designates the residual orbital-dependent interactions. In 
a real material, we expect V to vanish more quickly with distance than orbital-independent 
J. The charge ordering energy Ec[e\ = Ec[e\ is function of charge configuration e alone. Eq. ^ 
allows for consistent treatment of charge and orbital interactions. As will be discussed in our 
example, the separation into charge and orbital contributions is not unique and depends on 
our choice of independent parameters. If chosen appropriately, the parameters J and V can 
provide useful physical insight. They combine to derive E[e\, which is always meaningful. 
We will also see that lattice symmetry can further simplify eq. [2J 

B. Magnetite Structure 

We will concentrate on magnetite Fe304, a mixed- valence oxide with nominal iron valence 
between 2+ and 3+. At room temperature Fe304 has inverted cubic spinel structure FdSm 
with tetrahedral A sites occupied by one-third of the cations as Fe^^, and octahedral B 
sites by two-thirds of the cations with nominal valence 2.5-I-. At Ty ~ 120 K it undergoes 
the Verwey transition lowering the symmetrjAlES. Although the very existence of B-site 
Fe^"'"/Fe^"'" ordering at low-temperature (low-T) is not completely agreed upon, with some 
experimental data arguing against iiJEIll g^^id some supporting itP"^, we note that recent 
theoreticaP^"^ and experimentaP^HSS] results advocate both charge and orbital ordering in 
Fe304. The low-T structure has also been studied with model HamiltoniarP^'^. Piekarz 
and co-workers discussed the interplay of the electronic and ionic degrees of freedom in 
explaining the mechanism of the Verwey transition with first-principles and group theoretical 



arguments^Sl. 

Fig. [Ik shows the local environment of the Fe(B) atoms, where bonded Fe(B) and O 
atoms form corner-sharing cubes. The corner-sharing FeOg octahedra align as chains in the 
(110) directions (we always refer to the fee cell coordinates). Two parallel chains are shown 
in fig. ^ with nearest-neighbor (NN) and third-NN (3NN) B-sites highlighted. There are 
two distinct 3NN pairs: 3NNa on one chain with an intermediate Fe atom and 3NNb across 
two chains with no middle atom. The 2NN Fe(B) atoms do not share a {001} plane and are 
not shown. The B-sites form a pyrochlore lattice. Anderson found that the frustrated NN 






f*. 2 X(7it.-...3 ::Q 



7 




V 



FIG. 1. (Color online) a) Fragment of the Fe304 structure with Fe(B)-0 cubes, and the pyrochlore 
lattice formed by Fe(B) tetrahedra. b) One {001} plane of Fe atoms in (110) chains. Three distinct 
B pairs, two within a chain and one between two neighboring chains, are highlighted, c) Seven 
distinct t2g orbital interactions between coplanar Fe^+ or Fe'^"'' ions. Bond length difference in Fe^^ 
is exaggerated. Sphere indicates Fe^"*". 

interaction on this lattice leads to highly degenerate ground states with each B-tetrahedron 
occupied by two Fe^~^ and two Fe^+'^. The low-T CO structure proposed by Verwey^^'^ 
(for a picture see e.g. fig. 1 of ref. 123]) satisfies the Anderson condition, while some recent 
CO models, e.g. in refs. ITS f [T6 | andlSTj do not. To our understanding the low-temperature 
structure is still not fully resolvecP^'^. The low-T structure and the charge (and orbital) 
energetics therefore invites quantitative study. Here we present a detailed study of COO in 
Fe304 for the dual purpose of exploring first-principles calculation of orbital interactions in 
general, and to try to better understand the structure and origin of low-T phase. 



Magnetite is ferrimagnetic {T^ ~ 860 K) with antiparallel magnetic moments on the A 



and B sites at low-T. We fix the magnetic configuration as such in this work and focus on the 
charge and orbital degrees of freedom. The FeOg crystal field splits the five minority spin 
(i-states into three ^23 states and two higher-energy Cg states. At low temperature four states 
are accessible at each B site: the Fe^"*" hole and the Fe^~*" with t2g orbitals xy, yz, and xz (see 
fig.Ilfc). The symmetry and three-fold t2g degeneracy reduces the number of independent ^'s. 
We show in fig. ^ symmetrically distinct elements of the orbital interaction matrix V for 
NN (or 3NNa,b) pairs. Eq. |2]is simplified as follows: The electron chemical potential term Jj 
is unnecessary as the number of electrons is fixed in stoichiometric Fe304. The orbital point 
energy Vi is dropped due to the ^23 degeneracy. Some of the V matrix elements are linearly 
dependent on the J terms and may be removed. For example, V{1), orbital interaction 
type 1 (hole-hole) is already represented by the charge interaction Jiiij at the same sites. 
Only three matrix elements in fig. ^ are linearly independent within eq. [2j We choose to 
keep V{4), V{5) and V{7) and take V{6) as reference. The orbital-independent J is then 
unambiguously defined as between the reference orbital states, and V is the adjustment to J 
when the electronic states are not the reference. For example, the total interaction between 
NN xy electrons on the a6-plane is Jnat + VNAf(4). 

C. Computational details 

To parameterize the simplified eq. [2] we have performed Generalized Gradient Approxi- 
mation (GGA) + Hubbard U (GGA+U)^ calculations at f/cff = f/ — J = 4 eV (unless oth- 
erwise stated). All calculations were carried out using the VASP packag^^'^ with projector 
augmented wave (PAW) potential^, energy cutoff of 450 eV, and without any symmetry 
constraint on ionic and lattice relaxation. Each calculation was initialized in a specific con- 
figuration of charge and orbital order and self-consistently converged. We use supercells of 

TS^TS^-'-'TI^TI^^'-'-^-'-^-'-'-'-^-''^^ (designated I, II, III, IV) relative to the fee 
cell (see fig. 2 of ref. [16], where they were named P2/m, Fd3m, P2/c and Cc, respectively) 
and 2x1x1, 4x1x1 relative to the fee primitive cell. 

The issue or orbital moment is of considerable interest in the electronic structure of Fe304. 
Despite earlier reports of considerable orbital moment at the B-site Fe^"*" ions^, more recent 
measurements have found a relatively small orbital/spin moment raticPSESl. In this work we 
ignore spin-orbit coupling and assume completely suppressed orbital moment. 



The parameters in eq. |2] were determined with a iterative procedure commonly used 
in parametrization of cluster expansion: 1) fit GGA+U energy to eq. |2l 2) search with 
the obtained parameters for low-energy configurations, and 3) calculate new structures, if 
any, with GGA+U and go to step 1). This procedure was repeated until the parameters 
converged and no new ground state emerged. In the end, we calculated 365 distinct COO 
arrangements. 

III. RESULTS AND DISCUSSIONS 

In agreement with refs. [2nH231 and HOI charge disproportionation of ^ 0.2e is observed 
between Fe^"*" and Fe'^"'", with the t2g occupancies in the form of {a, 6, 6} or {6,6,6} {a ~ 
0.7 — 1.0, 6 ~ 0.0 — 0.3), respectively, clearly validating the notion of separating the Fe 
ions into distinct valence and orbital states. Another way to distinguish the ions is via the 
relaxed Fe-0 bond length. Typically the Fe^"'"-0 bond is 2.03 ±0.04 A, while the six Fe^"'"-0 
bonds are 2.11 A on average, with four elongated bonds of ~ 2.15 A within, and two shorter 
bonds perpendicular to, the plane of occupied orbital (see fig. [ip), proving that Fe^"*" is a 
Jahn- Teller active ion in Fe304. Our result confirms previous assessmentP^'^ that the Fe^"*" 
can be understood as Jahn- Teller active small polaron. 

The effect of orbital order on crystal structure is clearly demonstrated in fig. |2| where 
the lattice parameters (e.g. c) of 128 relaxed structures in the 1x1x1 fee cell are shown 
as function of /, the fraction of "perpendicular" (e.g. xy as opposed to c) orbitals among 
all t2g orbitals of B-site Fe^"*". Therefore a random configuration corresponds to / = 1/3. A 
linear fitting (dashed line in fig. u\ of (8.521 — 0.148/) A well captures the overall trend of 
the lattice parameters. The variation ~ 0.15 A across the range of / is considerably larger 
than the variation ~ 0.04 A at fixed /. 

The best fit of the GGA+U energies to eq. |2] has a cross validation score^ of 6 meV and 
root-mean-square (RMS) error of 5 meV per formula unit (f.u.), with 27 parameters. The 
orbital-independent J's include the constant term, 3 small triplet and 2 small quadruplet 
terms, and most significantly, 13 pair interactions shown in fig. [3j Note that these are ef- 
fective interactions including the many physical effects: electrostatics, screening, relaxation, 
covalency, etc. The NN pair ECI (solid line) is the largest orbital independent interaction 
(34 nieV), reflecting strong electrostatic repulsion. The orbital independent J's weaken con- 




0.25 0.5 0.75 

orbital fraction 



FIG. 2. (Color online) Lattice parameter versus the fraction of "perpendicular" orbitals. The gray 
dashed line is the best linear fit. 

siderably with distance and fall below 1 meV at > 10 A. A similar trend was observed in 
Li^FePOi^. 







C/eff = 4 eV 






UeS 


= 3.5 eV 




Pair 


y(4) 


V{5) 


y(7) 


F(4) 




V{5) 


^(7) 


NN 


106 


30 


3 


125 




33 


3 


3NNa 


32 


10 


-10 


38 




8 


-5 


3NNb 


6 








7 











TABLE I. Orbital interaction parameters V (fig. ITp) in meV at three Fe(B) pairs in fig. ITp. 

The orbital interactions V{n) [n = 4, 5, 7) for NN and 3NNa,b are listed in tablejl} VNAf(4) 
is by far the largest, which can be understood from 1) the a-bond-like orientation resulting 
in strong antiferromagnetic exchange interaction between same-spin electrons and 2) unfa- 
vorable quadrupole interactions. The orbital interactions obey the Kugel-Khomskii modeP, 
with V{6) = (reference) being the most stable. The 3NNa interactions are in general 
weaker than the NN, while the 3NNb and 2NN (not shown) are even smaller than 3NNa, 
though the distance is the same or shorter. Considering the different topology, the weak 

8 



- U-J=4eV 

- U-J=3.5eV 




6 8 

Pair distance/A 



FIG. 3. (Color online) Orbital-independent pair ECI versus pair distance at Ues = 4 eV (solid 
line) and 3.5 eV (dashed line). 

yet appreciable 3NNa orbital couplings suggest that the (110) chains of Fe/0 atoms may 
transmit exchange beyond NN. Note that a full interpretation of these effective interactions 
should include not only electronic but also lattice effects, e.g. Jahn- Teller coupling. 
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TABLE II. Calculated energy and predicted optimal orbital energy E^ (eq. pi) in meV/f.u. of 
certain structures. The first two structures are calculated with 00 optimized in different supercells 
(designated relative to fee), and the last two with COO patterns from ref. [23l The last column 
indicates whether the Anderson condition is met. 



To facilitate discussions we define the optimal orbital energy E^ of a given charge pattern 




FIG. 4. (Color online) The ground state COO pattern found in this work within supercell I 
(4= X 4= X 1 relative to fee) 

e minimized over OO's compatible with e. 



^™[?]=minE4e-l. 



(3) 



eGe 



The optimal orbital pattern is then defined as the one minimizing eq. |3} The Anderson 
degeneracy of the charge part of the energy landscape Ec, present with only NN interactions, 
is lifted by longer range charge interactions. The configurational space size for N f.u. is Cf^ . 
On the other hand, the search for the optimal orbital energy -E™[e] of a given CO e is also 
frustrated in a space of 3^. The orbital energy -^^[e] is also complicated by NN and longer 
range orbital interactions. Given the complex energy landscape, we use the above parameters 
to search the COO configuration space by enumeration in supercells II and III, and with 
Monte Carlo-based methods in larger supercells . The ground state COO pattern (fig. H| 
we find has the periodicity of 4= x 4= x 1 (though the periodicity of the structure is larger; 
see later discussions). As shown in fig. El the structure has equal number of Fe^"'"/Fe^"'" on 
each ab plane and uniform xz or yz electrons on alternate planes, i.e. no charge but orbital 
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modulation along c. We list the GGA+U energy of four structures, including our ground 
state, in table [Tij For the first two structures, the 00 was optimized in supercells I-IV to 
study their periodicity. First, consider the original Verwey CO model with alternate ab 
planes of electrons and holes, i.e. charge modulation along c. A large enough supercell is 
needed to find the optimal 00 of the Verwey CO model. In supercell I, all electrons occupy 
the xy state , while in larger cells II-IV, they equally occupy xy and xz to lower the energy 
by 9~11 meV (the variation is due to small convergence error in different supercells). Our 
ground state structure (fig. Ill) is confirmed to have the lowest GGA+U energy among all of 
our calculations. It has the same optimal 00 in the four supercells, but lowers its energy by 
15 meV in supercell III or IV compared to I or II, a difference too large to be a convergence 
error. It is found that in supercells I and II the Fe^"'"-0 bond lengths have similar distribution 
with standard deviation of 0.047 A, while in cells III and IV the deviation is 0.063 A. We 
believe that the energy difference has to do with lattice coupling of Jahn- Teller active Fe^"*": 
even with the same COO configuration the ionic positions in a small supercell are more 
constrained, reducing the chances of cooperative distortions. The structure therefore has 
the periodicity of supercell III, with space group P4i. Lastly, two previously proposed COO 
candidates^, with space group C(P^ and P2/c^^, respectively, are included. Both are less 
stable than our ground state. We confirm the conclusion of Jeng et al. that Cc is more 
stable than P2cp3. For the charge pattern of the Cc structure, an 00 2 meV lower than 
the one reported in ref. |231 is found. The predicted optimal orbital energy E^^ (eq. |3]) of 
our structure is 20 meV with all NN Fe^"'"-Fe^"'" interactions of the unfavorable type 5. Both 
P2/c and Cc have relatively small E^ of about 10 meV, since many of their NN interactions 
are of type 6 or 7. 

The experimental low-temperature structure of magnetite is not yet known clearly enough 
to compare with, but the fact that our ground state has different unit cell and space group 
than found in experiment!^ means we have not completely resolved the problem. Neverthe- 
less, several observations can be made from our results. First, like the charge energy, the 
orbital energy is also frustrated. For example, the optimal 00 for each of the four structures 
in table [n] include unfavorable orbital interactions VNNin) {n = 5, 7). Secondly, the charge 
and orbital energies are competing: a structure may have low charge energy E^ or orbital 
energy E^, but not both. Taking the Anderson condition as an approximate indicator of low 
Ec, our ground state COO has low E^. but relatively large E'^, while the Cc and P2/c are the 
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opposite (table [n]). Thirdly, the frustrated, competing interactions make the ground state 
search sensitive to the interaction parameters. It is possible that our search failed to find 
a ground state in a larger supercell because of the numerical sensitivity. Other possibilities 
for the incomplete agreement with experimental supercell might be the missing physics not 
described by our calculations. This includes spin fluctuations beyond the assumption of a 
fixed magnetic configuration, and may also include fluctuations between complex structures 
of close energies. Lastly, our results illustrate some possible mechanisms that can break the 
cubic symmetry and form the low-T GS structure: (1) charge order as Verwey originally 
proposed, (2) charge and orbital order as exemplified in the Verwey CO model whose COO 
supercell is larger than the CO supercell, and (3) lattice coupling of Fe^^ ions as seen in 
our structure (fig. Ill), the periodicity of which, decided by the arrangement of Jahn- Teller 
distortions, is larger than that of the COO. 

To evaluate the impact of the Hubbard term Ues in our first-principles approach, we have 
calculated 300 structures with a smaller Ucs = 3.5 eV. The pair ECIs for Ucs = 3.5 eV 
are shown as the dashed line in fig. |3} They are slightly reduced compared to Ues = 4 eV, 
mainly because the largely electrostatic ECIs scale as (Ag)^ where Aq is the charge difference 
between the 2+/3+ ions. With smaller Ues the d-electrons become more delocalized and 
Aq generally decreaseP^^ll. As shown in table |l] the orbital interactions V are relatively 
stable yet some are notably larger at Ucs = 3.5 eV. Presumably the reason is that the 
exchange integral, sensitive to the spatial distribution of the wavef unctions, increases with 
delocalization. It is also possible that the mathematical separation in eq. [2] of charge and 
orbital terms is not physically complete, and there is some compensation in J and V with 
varying Ues- The smaller Hubbard parameter does not considerably change the results in 
table HI] and related discussions. 

IV. CONCLUSIONS 

In this work we have attempted to describe the charge and orbital degrees of freedom in 
Fe304 with a classical effective energy model. Electronic and lattice effects are both included 
through first-principles calculated energies from which the model is parametrized. The 
calculated charge and orbital interactions in Fe304 are found to be physically meaningful. 
The energy landscape is complex in terms of frustrated charge and orbital interactions 
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as well as their competition. Additionally, although our predicted ground state structure 
has smaller periodicity than experimentally observed, it reveals the possibility that not 
only charge and orbital ordering, but the Jahn- Teller lattice distortions may also decide 
the structure. Therefore this work may help better understand the problem of the low-T 
magnetite structure. Beyond magnetite, our approach can be easily adapted to explore other 
transition metal oxides where charge and/or orbital order exist. 

This work is supported by DOE under contract DE-FG02-96ER45571, and in part by 
NSF through the National Partnership for Advanced Computing Infrastructure using SDSC 
Datastar. FZ thanks Dr. C Fischer for his help in data analysis. 
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